/* Colluding Against Workers: Main model: bootstrapping
	Delabastita & Rubens 
------------------------------------------*/

clear
set more off
cd "$station"
   
   
* global bs = 2 
local B = $bs		 // number of bootstrap iterations (set this to 200)
local C = 50			// convergence max iterations (set this to 50)
 
set matsize 1000	
set more off
clear  

** Matrices to store estimates
 
forvalues m = 1(1)3 {
foreach mod in   "bb" "ols" {
foreach coef in "thl_`mod'" "thm_`mod'"  "scale_`mod'" "bl_`mod'" "rho_`mod'"  "bk_`mod'" "blk_`mod'" "bmk_`mod'"   "blm_`mod'"  "bm_`mod'" "bl2_`mod'" "bm2_`mod'" "bk2_`mod'" {
matrix _`coef'`m' = J(`B', 1, 0)
matrix m_`coef'`m' = J(`B', 1, 0)
}
foreach coef in "md_bb" "mu_bb"  "conv"    {
matrix _`coef'`m' = J(`B', 1, 0)
matrix m_`coef'`m' = J(`B', 1, 0)
matrix a_`coef'`m' = J(`B', 1, 0)
}
}
}

* Matrix to store returns to scale model

foreach coef in "bl" "bm" "bk" "md" "mu"  {
forvalues n = 0(5)15 {
	matrix _`coef'_nu`n' = J(`B', 1, 0)
	matrix m_`coef'_nu`n' = J(`B', 1, 0)
}	
}

* Matrix to store pf coefs by time blocks

foreach coef in "bl" "bm" "bk" {
forvalues n = 1 / 2 {
	matrix _`coef'_blockbis`n' = J(`B', 1, 0)
}	
}

forvalues m = 1(1)4 {
matrix _kappa_dq`m' = J(`B', 1, 0)
matrix m_kappa_dq`m' = J(`B', 1, 0)
}

* IV comparisons 
foreach coef in "bl" "bm" "bk" {
forvalues n = 1 / 3 {
	matrix _`coef'_iv`n' = J(`B', 1, 0)
}	
}


* Matrix to store rev pf

foreach coef in "mu" "md"  {
forvalues n = 1 / 12 {
	matrix _`coef'_bb`n' = J(`B', 1, 0)
	matrix m_`coef'_bb`n' = J(`B', 1, 0)
	matrix a_`coef'_bb`n' = J(`B', 1, 0)
}	
}


foreach coef in "bl" "bm" "bk" "rho" "rho1" "rho2" {
forvalues n = 1 / 12 {
	matrix _`coef'_bb`n' = J(`B', 1, 0)
}	
}

* Matrix to store markup response

forvalues n = 1 / 2 {
	matrix _mucar`n' = J(`B', 1, 0)
}	


foreach coef in "blt_bb5" "bl_bb5" "bt_bb5" "bnu" {
matrix _`coef' = J(`B', 1, 0)	
}

* Matrix to store CS-based coefficients

foreach coef in "bl" "bm" "bk" "mu" {
matrix _`coef'_cs_noc = J(`B', 1, 0)
matrix _`coef'_cs_col = J(`B', 1, 0)
}	

  
 
forvalues m = 1/2 {
foreach coef in   "th_rail`m'" "th_union`m'" "th_tram`m'"   "th_nm1`m'" "th_nm2`m'" "th_nm3`m'" "th_urban`m'"  "th_car`m'"   ///
 "th_lib`m'" "th_soc`m'" "th_com`m'" "th_oth`m'" "th_libp`m'" "th_socp`m'" "th_comp`m'" "th_othp`m'"   { 
matrix _`coef' = J(`B', 1, 0)
matrix m_`coef' = J(`B', 1, 0)
} 
}

  
use ./data/temp/tempap, clear 
 di _N 
 
** Start bootstrapping iteration


forvalues b = 1 /`B' {
       preserve
         set seed `b'
          bsample   ,  cluster(mineid) idcluster(idb)
         drop mineid
           rename idb mineid   

  /* A. Production function: robustness checks */

  ** 1/ Baseline model

drop *bb1

gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm  L.lk lk L.lsegers  )  conv_maxiter(`C') 
 
mat coef = e(b)
mat list coef
gen rho_bb1 = coef[1,1]
gen bl_bb1 = coef[1,2]
gen bm_bb1 = coef[1,3]
gen bk_bb1 = coef[1,4]
gen bc_bb1 = coef[1,5]
gen blm_bb1  = .
gen blk_bb1 = .

gen thl_bb1 = bl_bb1  
gen thm_bb1 = bm_bb1  

gen scale_bb1 = thl_bb1 +  thm_bb1 + bk_bb1 //+ bd_bb1 


** 1b. Baseline model fixing returns to scale 

forvalues n = 0(5)15 {
local nu = 1 +`n'/100
gen nu`n' = `nu'
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -(`nu'-{bl}-{bm})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm  L.lk lk L.lsegers  )  conv_maxiter(50) 
 
mat coef = e(b)
mat list coef
gen rho_nu`n' = coef[1,1]
gen bl_nu`n' = coef[1,2]
gen bm_nu`n' = coef[1,3]
gen bk_nu`n' = 1-bl_nu`n'-bm_nu`n'
gen bc_nu`n' = coef[1,4]
gen blm_nu`n'  = .
gen blk_nu`n' = .

gen thl_nu`n' = bl_nu`n' 
gen thm_nu`n' = bm_nu`n'  

gen mu_nu`n' = bm_nu`n' /am 
gen md_nu`n' = (bl_nu`n' /al)/mu_nu`n' 

}
 

** 2/ Interaction effect materials and labor

gen lemplk = lemp*lk
gen lwmlk = lwm*lk

* OLS

reg lq lemp lwm  lk lemplk  lwmlk

mat coef = e(b)
mat list coef
gen bl_ols2 = coef[1,1]
gen bm_ols2 = coef[1,2]
gen bk_ols2 = coef[1,3]
gen blk_ols2 = coef[1,4]
gen bmk_ols2 = coef[1,5]
gen bc_ols2 = coef[1,6]

gen thl_ols2 = bl_ols2   + blk_ols2 * lk
gen thm_ols2 = bm_ols2   
gen scale_ols2 = thl_ols2 +  thm_ols2 + bk_ols2 

local r2 = string(e(r2)) 		 
local r2_ols2 = substr("`r2'" ,1,4) 
 
* GMM

gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk) ///
 -({blk})*(lemplk - {rho}*L.lemplk)   -({bmk})*(lwmlk - {rho}*L.lwmlk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})    )      , ///
inst(  L.lemp  L.lwm L.lemplk L.lwmlk  lsegers L.lk lk  )  conv_maxiter(`C') 


mat coef = e(b)
mat list coef
gen rho_bb2 = coef[1,1]
gen bl_bb2 = coef[1,2]
gen bm_bb2 = coef[1,3]
gen bk_bb2 = coef[1,4]
gen blk_bb2 = coef[1,5]
gen bmk_bb2 = coef[1,6]
gen bc_bb2 = coef[1,7]

gen thl_bb2 = bl_bb2  + blk_bb2 * lk
gen thm_bb2 = bm_bb2   + bmk_bb2 * lk
gen scale_bb2 = thl_bb2 +  thm_bb2 + bk_bb2  

 
gen ltfp2 = lq - bl_bb2*lemp-bk_bb2*lk - bc_bb2 -bm_bb2*lwm  - blk_bb2*lemplk  
gen tfp2 = exp(ltfp2)
   
** 3/ Translog production function  

xtset mineid yr
gen Llklemp =  lk*L.lemp
gen Llklwm =  lk*L.lwm
 gen lemplwm = lemp*lwm
foreach var of varlist lemp lwm lk {
 gen `var'2 = `var'^2	
}
   
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)     -({bk})*(lk - {rho}*L.lk) ///
  -({blk})*(lemplk - {rho}*L.lemplk)  -({blm})*(lemplwm - {rho}*L.lemplwm) -({bmk})*(lwmlk - {rho}*L.lwmlk) ///		 
 -({bl2})*(lemp2 - {rho}*L.lemp2) -({bm2})*(lwm2 - {rho}*L.lwm2) -({bk2})*(lk2 - {rho}*L.lk2)     ///		 
   -{bc}*(1 - {rho})         ) , ///
inst(  L.lemp  L.lwm L.lk lk  L.lsegers  L.lemplk L.lemplwm L.lwmlk   Llklemp Llklwm  L.lemp2 L.lwm2 lk2 L.lk2)  conv_maxiter(`C') 
 
 
 
mat coef = e(b)
mat list coef

gen rho_bb3 = coef[1,1]  
gen bl_bb3 = coef[1,2]
gen bm_bb3 = coef[1,3]
gen bk_bb3 = coef[1,4]
gen blk_bb3 = coef[1,5]
gen blm_bb3 = coef[1,6]
gen bmk_bb3 = coef[1,7]
gen bl2_bb3 = coef[1,8]
gen bm2_bb3 = coef[1,9]
gen bk2_bb3 = coef[1,10]

gen thl_bb3 = bl_bb3 + blm_bb3 * lwm + blk_bb3 * lk + bl2_bb3 * 2 *lemp 
gen thm_bb3 = bm_bb3  + blm_bb3 * lemp + bmk_bb3*lk + bm2_bb3 * 2*lwm 
  

** 4/ Time-varying production function

* Two time blocks (1845-1879, 1880-1913)	 

gen blockbis1 = yr >=1845 & yr<=1879
gen blockbis2 = yr >=1880 & yr<=1913
 
gen bl_bb4 = .
gen bm_bb4 = .
gen bk_bb4 = .
gen rho_bb4 = .

 
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)       ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})    )    if blockbis1==1  , ///
inst(  L.lemp  L.lwm   L.lsegers   L.lk lk  )  conv_maxiter(`C') //technique(bfgs)	// does converge using BFGS

mat coef = e(b)
mat list coef
replace rho_bb4 = coef[1,1] if blockbis1==1
replace bl_bb4 = coef[1,2] if blockbis1==1
replace bm_bb4 = coef[1,3] if blockbis1==1 
replace bk_bb4 = coef[1,4] if blockbis1==1

gen bl_blockbis1 = coef[1,2]
gen bm_blockbis1 = coef[1,3]
gen bk_blockbis1 = coef[1,4]
 
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)       ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})    )    if blockbis2==1  , ///
inst(  L.lemp  L.lwm    L.lsegers  L.lk lk  )  conv_maxiter(`C') //technique(bfgs)

mat coef = e(b)
mat list coef
replace rho_bb4 = coef[1,1] if blockbis2==1
replace bl_bb4 = coef[1,2] if blockbis2==1
replace bm_bb4 = coef[1,3] if blockbis2==1
replace bk_bb4 = coef[1,4] if blockbis2==1

gen bl_blockbis2 = coef[1,2]
gen bm_blockbis2 = coef[1,3]
gen bk_blockbis2 = coef[1,4]
 
 
gen thl_bb4 = bl_bb4  
gen thm_bb4 = bm_bb4  

 
** 5/ Time-varying coefficients: interaction effects on beta_l

gen lempyr= lemp*yr 
gen lwmyr = lwm*yr 
gen lkyr = lk*yr
   
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)       ///		// equation (8) with blk = 0
 -{blt}*(lempyr - {rho}*L.lempyr)         -{bt}*(yr - {rho}*L.yr)      -{bc}*(1 - {rho})    )     , ///
inst(  L.lemp  L.lwm   L.lsegers    L.lk lk yr L.lempyr     )  conv_maxiter(`C')  //technique(bfgs)
mat coef = e(b)
mat list coef
gen rho_bb5 = coef[1,1]  
gen bl_bb5 = coef[1,2]  
gen bm_bb5 = coef[1,3]   
gen bk_bb5 = coef[1,4]   
gen blt_bb5 = coef[1,5] 
gen bt_bb5 = coef[1,6] 
 
gen thl_bb5 = bl_bb5 + blt_bb5*yr
gen thm_bb5 = bm_bb5  
 

** 6/ Add linear time trend
 
gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)   -({bt})*(yr - {rho}*L.yr)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})    )      , ///
inst(  L.lemp  L.lwm  L.lsegers    L.lk lk  yr)  conv_maxiter(`C') //technique(bfgs)

mat coef = e(b)
mat list coef

gen rho_bb6 = coef[1,1]  
gen bl_bb6 = coef[1,2]  
gen bm_bb6 = coef[1,3]   
gen bk_bb6 = coef[1,4]  

gen thl_bb6 = bl_bb6  
gen thm_bb6 = bm_bb6  


** 7/ Unobserved input prices

gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)       ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})        -{bp}*(lp - {rho}*L.lp) )     , ///
inst(  L.lemp  L.lwm    L.lsegers     L.lk lk  L.lp lp)  conv_maxiter(`C') //technique(bfgs)
    
mat coef = e(b)
mat list coef

gen rho_bb7 = coef[1,1]  
gen bl_bb7 = coef[1,2]  
gen bm_bb7 = coef[1,3] 
gen bk_bb7 = coef[1,4] 
  
gen thl_bb7 = bl_bb7  
gen thm_bb7 = bm_bb7 

** 8/ Revenue production function

gen lrev= log(p*q)

gmm (lrev - {rho}*L.lrev -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)       ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})         )     , ///
inst(  L.lemp  L.lwm   L.lsegers    L.lk lk  )  conv_maxiter(`C') //technique(bfgs)
   
mat coef = e(b)
mat list coef
gen rho_bb8 = coef[1,1]  
gen bl_bb8 = coef[1,2]  
gen bm_bb8 = coef[1,3] 
gen bk_bb8 = coef[1,4] 
   
gen thl_bb8 = bl_bb8  
gen thm_bb8 = bm_bb8 
   
   
** 9/ Quality control
 
foreach var of varlist q_leans_tot_t   {
gen s`var' = `var'/q
gen ls`var'= log(s`var')
}
replace sq_leans_tot_t = 0 if sq_leans_tot_t ==.

reg ltfp1 lsq_lean 
xtreg ltfp1 lsq_lean , fe r

gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)       ///		// equation (8) with blk = 0
 -{bs}*(sq_leans_tot_t - {rho}*L.sq_leans_tot_t)    -{bc}*(1 - {rho})         )     , ///
inst(  L.lemp  L.lwm    L.lsegers    L.lk lk sq_leans_tot_t L.sq_leans_tot_t )  conv_maxiter(`C') //technique(bfgs)
  
 
mat coef = e(b)
mat list coef
gen rho_bb9 = coef[1,1]  
gen bl_bb9 = coef[1,2]  
gen bm_bb9 = coef[1,3] 
gen bk_bb9 = coef[1,4] 
 
gen thl_bb9 = bl_bb9  
gen thm_bb9 = bm_bb9

 
*gen p98 = yr>=1897


** 10/ serial correlation = 1
  
gmm (lq -  L.lq -{bl}*(lemp - L.lemp)     -({bm})*(lwm -  L.lwm)   -({bk})*(lk -  L.lk)       ///		// equation (8) with blk = 0
           )     , ///
inst(  L.lemp  L.lwm  L.lsegers   L.lk lk   )  conv_maxiter(`C') // technique(bfgs)
  
mat coef = e(b)
mat list coef
  
gen rho_bb10 = 1
gen bl_bb10 = coef[1,1]  
gen bm_bb10 = coef[1,2] 
gen bk_bb10 = coef[1,3]   
 
gen thl_bb10 = bl_bb10  
gen thm_bb10 = bm_bb10


** 11/ Estimate production function without agricultural wage 

gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm  L.lk lk   )  conv_maxiter(`C') technique(bfgs)

mat coef = e(b)
mat list coef
gen rho_bb11 = coef[1,1]
gen bl_bb11 = coef[1,2]
gen bm_bb11 = coef[1,3]
gen bk_bb11 = coef[1,4]
gen bc_bb11 = coef[1,5]
gen blm_bb11  = .
gen blk_bb11 = .

gen thl_bb11 = bl_bb11  
gen thm_bb11 = bm_bb11  

 
gen scale_bb11 = thl_bb11 +  thm_bb11 + bk_bb11 
 

 * estimate with agricultural shocks themselves. 

gen limp_p_isic = log(imp_p_isic)	// log import price of agricultural products 
gen limp_p_keygrains = log(imp_p_keygrains)	// log import price of key grain products 
gen lprod_agri = log(agriculturalproduct/agr_emp_imputed)

* i. main IV

gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm  L.lk lk  L.lsegers  )  conv_maxiter(50)  
mat coef = e(b)
mat list coef
gen rho_iv1 = coef[1,1]
gen bl_iv1 = coef[1,2]
gen bm_iv1 = coef[1,3]
gen bk_iv1 = coef[1,4]

* ii. no instrument 

gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm  L.lk lk   )  conv_maxiter(50)   technique(bfgs)
mat coef = e(b)
mat list coef
gen rho_iv2 = coef[1,1]
gen bl_iv2 = coef[1,2]
gen bm_iv2 = coef[1,3]
gen bk_iv2 = coef[1,4]


* iii. agricultural import prices

gmm (lq - {rho}*L.lq -{bl}*(lemp - {rho}*L.lemp)     -({bm})*(lwm - {rho}*L.lwm)   -({bk})*(lk - {rho}*L.lk)      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho})      )      , ///
inst(  L.lemp  L.lwm  L.lk lk   L.limp_p_keygrains   L.lprod_agri )  conv_maxiter(50)  

mat coef = e(b)
mat list coef
gen rho_iv3 = coef[1,1]
gen bl_iv3 = coef[1,2]
gen bm_iv3 = coef[1,3]
gen bk_iv3 = coef[1,4]


 
 
 
  
** 12/ AR(2) production function

gmm (lq - {rho}*L.lq - {rho2}*L2.lq -{bl}*(lemp - {rho}*L.lemp - {rho2}*L2.lemp)     -({bm})*(lwm - {rho}*L.lwm - {rho2}*L2.lwm)   -({bk})*(lk - {rho}*L.lk - {rho2}*L2.lk) ///
      ///		// equation (8) with blk = 0
   -{bc}*(1 - {rho} - {rho2})    )      , ///
inst(  L.lemp L2.lemp  L.lwm L2.lwm    L.lsegers L2.lsegers L.lk lk  L2.lk )  conv_maxiter(`C') //technique(bfgs)

mat coef = e(b)
mat list coef
drop *bb12

gen rho1_bb12 = coef[1,1]
gen rho2_bb12 = coef[1,2]
gen bl_bb12 = coef[1,3]
gen bm_bb12 = coef[1,4]
gen bk_bb12 = coef[1,5]
gen bc_bb12 = coef[1,6]
 
gen thl_bb12 = bl_bb12  
gen thm_bb12 = bm_bb12  

 
** Compute markdowns and markups 
   
  forvalues m = 1/12 {
gen mu_bb`m' = thm_bb`m' /am 
gen md_bb`m' = (thl_bb`m' /al)/mu_bb`m' 
}

drop md_bb8 mu_bb8 	// rev pf - markdowns computed differently and markup not identified

gen mdl_bb8 = (bl_bb8 /al) 
gen md_bb8 = mdl_bb8
gen mu_bb8 = .
gen mdm_bb8 = (bm_bb8 /am)

 
** Other robustness checks

* Productivity innovations serially correlated?

xtset mineid yr
gen ltfp_bb1 = lq - bl_bb1*lemp-bk_bb1*lk - bc_bb1 -bm_bb1*lwm  

drop nu 
gen nu =  ltfp_bb1 - rho_bb1*L.ltfp_bb1 

reg nu L.nu
gen bnu = _b[L.nu]

* What happens if \rho = 1?

gmm (lq -  L.lq -{bl}*(lemp -  L.lemp)     -({bm})*(lwm -L.lwm)   -({bk})*(lk -  L.lk)       ///		// equation (8) with blk = 0
   )      , ///
inst(  L.lemp  L.lwm   L.lsegers      L.lk lk  )  conv_maxiter(`C') // technique(bfgs)

 
 
* interaction effects

foreach mod in "ols" "bb" {
foreach var in "bl" "bm" "bk"  "blk"   "bmk" "thm" "thl" "scale"{
gen `var' = `var'_`mod'2
}
estpost su bl bm bk     blk bmk
est store pf_`mod'2
estpost su thl thm scale 
est store th_`mod'2
drop bl bm bk    thl scale thm   blk bmk
}

  

* Cost shares approach - robustness check
 
gen rts = 1	// returns to scale parameter
gen bl_cs_col = rts*(wl*psim_dq4/(wm + wl*psim_dq4 + winv)) 
gen bm_cs_col = rts*(wm/(wm + wl*psim_dq4 + winv)) 
gen bk_cs_col = rts - bl_cs_col - bm_cs_col

gen bl_cs_noc = rts*(wl*psic_dq4/(wm + wl*psic_dq4+ winv)) 
gen bm_cs_noc = rts*(wm/(wm + wl*psic_dq4+ winv)) 
gen bk_cs_noc = rts - bl_cs_noc - bm_cs_noc
 
sum bl_cs_col bl_cs_noc, d
sum bm_cs_col bm_cs_noc, d

gen mu_cs_noc = bm_cs_noc / am
gen mu_cs_col = bm_cs_col / am

sum mu_cs*
  
sum pop
bys dq4id yr: egen pop_dq4 = mean(pop) 	// pop comes at dq4 level
drop emp_dq4
bys dq4id yr: egen nemp_dq4 = sum(emp/ndays)

gen nemp = emp/ndays

drop semp_dq4 lsemp_dq4
gen semp_dq4 = nemp/pop_dq4
gen lsemp_dq4 = log(nemp/pop_dq4)
gen oo_dq4 = (pop_dq4-nemp_dq4) /pop_dq4
gen loo_dq4  = log(oo_dq4)

gen lhs_dq4  = lsemp_dq4 - loo_dq4 
drop uw
gen uw = wl/emp		// daily wage 
gen uwn = wl/nemp	// annual wage
reg lhs_dq4   uw
drop luw
gen luw = log(uw)
gen luwn = log(uwn)

gen demshock = yr>= 1871 & yr<= 1875
gen p98 = yr>=1898
gen e_car_p98 = e_car*p98


ivregress 2sls  lhs_dq4 (uw   = demshock e_car_p98) e_car p98  , r
drop bw
gen bw = _b[uw]

gen md_lin = (bw*uw*(1-semp_dq4))^(-1)+1
sum md_lin, d

ivregress 2sls  lhs_dq4 (luw = demshock e_car_p98) e_car p98  , r
gen blw = _b[luw]
gen md_con = (blw*(1-semp_dq4))^(-1)+1
sum md_con, d

* Fully collusive markdown: need to solve for equilibrium wage under full collusion

gen kappa_dq4_lin = md_bb1 / md_lin
gen kappa_dq4_con = md_bb1 / md_con

sum kappa_dq4_lin kappa_dq4_con, d

foreach var of varlist kappa_dq4_con kappa_dq4_lin {
bys yr: egen med`var' = median(`var')	
bys yr: egen av`var' = mean(`var')	
}


/* C. Markups and markdowns
--------------------------------*/
  
  

forvalues m = 2(1)2 {
foreach mod in   "bb" "ols" {
foreach var in  "thm_`mod'" "thl_`mod'" "bk_`mod'" "bl_`mod'"   "bm_`mod'"  /*"bd_`mod'" */ "scale_`mod'" {
 egen av`var'`m' = mean(`var'`m')
}
}
}



forvalues m = 1(1)12 {
foreach var in  "md_bb" "mu_bb"   {
bys yr: egen  med`var'`m' = median(`var'`m')
bys yr: egen  av`var'`m' = mean(`var'`m')
gen s`var'`m' = `var'`m'*sempag
bys yr: egen  ag`var'`m' = sum(s`var'`m')
}
}

 
foreach var in "md" "mu"  {
forvalues n = 0(5)15 {
  egen  m`var'_nu`n' = median(`var'_nu`n')
}
}

sum mmd_nu5 , d 
 

 forvalues m = 12(1)12 {
foreach var in   "bk" "bl"   "bm"  "rho2" "rho1" {
 egen av`var'_bb`m' = mean(`var'_bb`m')
}
}


 
* Robustness: markups and cartel
gen lmu_bb1 = log(mu_bb1)
reg lmu_bb1 e_car e_car_p98 p98 	// yes markup increases at cartel firms after 1898
gen mucar1 = _b[e_car_p98] 
areg lmu_bb1 e_car e_car_p98 p98 , absorb(mine_id)	// yes markup increases at cartel firms after 1898
gen mucar2 = _b[e_car_p98] 

** Cost dynamics

xtset mineid yr
gen qcum = 0
forvalues m = 1/72 {
replace qcum = qcum+ L`m'.q if L`m'.q~=. 
}

gen lqcum = log(qcum)
label var lqcum "log(Cum. output)"
label var ltfp1 "log(TFP)"
 
bys mineid: egen startyr = min(yr)
gen age = yr-startyr
gen lage = log(age)
reg ltfp1 lqcum lage i.yr, r


keep in 1/`B'
 
 
forvalues m = 1(1)3 {
foreach var in    "md_bb"  "mu_bb"   {
 mkmat av`var'`m', matrix(av`var'`m')
  mkmat med`var'`m', matrix(med`var'`m')
  mkmat ag`var'`m', matrix(ag`var'`m')
 matrix _`var'`m'[`b',1] = av`var'`m'[1,1]			// store average markdown for this bootstrap iteration in matrix
  matrix m_`var'`m'[`b',1] = med`var'`m'[1,1]		 
  matrix a_`var'`m'[`b',1] = ag`var'`m'[1,1]		 
}  
}


* IV comparisons 
foreach coef in "bl" "bm" "bk" {
forvalues n = 1 / 3 {
 mkmat `coef'_iv`n', matrix(`coef'_iv`n')
matrix _`coef'_iv`n'[`b',1] = `coef'_iv`n'[1,1]
}	
}


forvalues m = 4(1)12 {
foreach var in    "md_bb"  "mu_bb"   {
 mkmat av`var'`m', matrix(av`var'`m')
 matrix _`var'`m'[`b',1] = av`var'`m'[1,1]			// store average markdown for this bootstrap iteration in matrix
}  
}

forvalues m = 1(1)2 {
 mkmat mucar`m', matrix(avmucar`m')
 matrix _mucar`m'[`b',1] = avmucar`m'[1,1]			// store average markdown for this bootstrap iteration in matrix 
}

 
foreach mod in   "bl" "bm" "bk" "rho"  {
forvalues m = 1(1)11 {
 mkmat `mod'_bb`m', matrix(`mod'_bb`m')
 matrix _`mod'_bb`m'[`b',1] = `mod'_bb`m'[1,1]			// store average markdown for this bootstrap iteration in matrix 
}
}

foreach mod in  "bl" "bm" "bk" "rho1"  "rho2"  {
forvalues m = 12(1)12 {
 mkmat `mod'_bb`m', matrix(`mod'_bb`m')
 matrix _`mod'_bb`m'[`b',1] = `mod'_bb`m'[1,1]			// store average markdown for this bootstrap iteration in matrix 
}
}
 
 
forvalues m = 1(1)2 {
foreach mod in   "bl" "bm" "bk" {
 egen av`mod'_blockbis`m' = mean(`mod'_blockbis`m')
 mkmat av`mod'_blockbis`m', matrix(av`mod'_blockbis`m')
 matrix _`mod'_blockbis`m'[`b',1] = av`mod'_blockbis`m'[1,1]		 
}
}

foreach mod in   "bt_bb5" "bl_bb5" "blt_bb5" "bnu" {
 egen av`mod' = mean(`mod')
 mkmat av`mod', matrix(av`mod')
 matrix _`mod'[`b',1] = av`mod'[1,1]		 
}
 

forvalues m = 3(1)3 {
foreach mod in   "bl" "bm" "bk" "blk" "bmk" "blm" "bk2" "bl2" "bm2" {
 egen av`mod'_bb`m' = mean(`mod'_bb`m')
 mkmat av`mod'_bb`m', matrix(av`mod'_bb`m')
 matrix _`mod'_bb`m'[`b',1] = av`mod'_bb`m'[1,1]		 
}
}

foreach coef in "bl" "bm" "bk"  {
forvalues n = 0(5)15 {
egen av`coef'_nu`n' = mean(`coef'_nu`n')
mkmat av`coef'_nu`n', matrix(av`coef'_nu`n')
matrix _`coef'_nu`n'[`b',1] = av`coef'_nu`n'[1,1]		 
}	
}
foreach coef in "md" "mu"  {
forvalues n = 0(5)15 {
egen av`coef'_nu`n' = mean(`coef'_nu`n')	// take median for md, mu 
mkmat av`coef'_nu`n', matrix(av`coef'_nu`n')
matrix _`coef'_nu`n'[`b',1] = av`coef'_nu`n'[1,1]		 
egen med`coef'_nu`n' = mean(m`coef'_nu`n')	// take median for md, mu 
mkmat med`coef'_nu`n', matrix(m`coef'_nu`n')
matrix m_`coef'_nu`n'[`b',1] = m`coef'_nu`n'[1,1]		 
}	
}



forvalues m = 1/4 {
 mkmat avkappa_dq`m', matrix(avkappa_dq`m')
 mkmat medkappa_dq`m', matrix(medkappa_dq`m')
 }
 
   
forvalues m = 2(1)2 {
foreach mod in   "bb" "ols" {
foreach var in "thl_`mod'"  "thm_`mod'" "bk_`mod'" "bl_`mod'"    "bm_`mod'" "scale_`mod'"   {
mkmat av`var'`m', matrix(av`var'`m')
matrix _`var'`m'[`b',1] = av`var'`m'[1,1]			// store average markdown for this bootstrap iteration in matrix
}
}
}
 
foreach mod in   "bl" "bm" "bk" "mu" {
egen av`mod'_col = mean(`mod'_cs_col)
egen av`mod'_noc = mean(`mod'_cs_noc)
mkmat av`mod'_col, matrix(av`mod'_col)
mkmat av`mod'_noc, matrix(av`mod'_noc)
matrix _`mod'_cs_col[`b',1] = av`mod'_col[1,1]		 
matrix _`mod'_cs_noc[`b',1] = av`mod'_noc[1,1]		 
}
 

 
 
restore			// end bootstrapping loop here
}

 
 
 
**********

forvalues m = 1(1)12 {
foreach var in "md_bb"   "mu_bb"    {
svmat _`var'`m'   
svmat m_`var'`m' 
svmat a_`var'`m' 
}
}


foreach coef in "bl" "bm" "bk" "md" "mu"{
forvalues n = 0(5)15 {
svmat _`coef'_nu`n'   
svmat m_`coef'_nu`n'   
}	
}


* IV comparisons 
foreach coef in "bl" "bm" "bk" {
forvalues n = 1 / 3 {
svmat _`coef'_iv`n' 
}	
}


forvalues m = 1(1)4 {
foreach var in "kappa"    {
svmat _`var'_dq`m'   
svmat m_`var'_dq`m'   
 }
}
 
 forvalues m = 1(1)2 {
svmat _mucar`m' 
}


forvalues m = 1(1)2 {
foreach mod in   "bb" "ols" {
foreach var in "rho_`mod'" "thl_`mod'"  "thm_`mod'" "bk_`mod'" "bl_`mod'" "bm_`mod'" "blk_`mod'" "bmk_`mod'"  /*"bd_`mod'"*/   "scale_`mod'" /*"blloc_`mod'" "blcut_`mod'"*/ {
svmat _`var'`m'   
}
}
}


forvalues m = 3(1)3 {
foreach mod in   "bl" "bm" "bk" "blk" "bmk" "blm" "bk2" "bl2" "bm2" {
svmat _`mod'_bb`m'   
}
}


foreach mod in   "bl" "bm" "bk" "rho"   {
forvalues m = 4(1)11 {
 svmat _`mod'_bb`m'  			// store average markdown for this bootstrap iteration in matrix 
}
}
 
foreach mod in   "bl" "bm" "bk" "rho1" "rho2"   {
forvalues m = 12(1)12 {
 svmat _`mod'_bb`m'  			// store average markdown for this bootstrap iteration in matrix 
}
}
 
  
 
 
foreach mod in   "bl" "bm" "bk"  {
forvalues m = 1(1)2 {
svmat _`mod'_blockbis`m'   
}
}
  
foreach mod in   "bl_bb5" "blt_bb5" "bt_bb5" "bnu"   {
capture: svmat _`mod' 
}

foreach mod in   "bl" "bm" "bk" "mu" {
svmat _`mod'_cs_noc   
svmat _`mod'_cs_col  
}
  
forvalues m = 1/2 {
foreach coef in   "th_rail`m'" "th_tram`m'"   "th_nm1`m'" "th_nm2`m'" "th_nm3`m'" "th_urban`m'" "th_car`m'"  ///
 "th_lib`m'" "th_soc`m'" "th_com`m'" "th_oth`m'" "th_libp`m'" "th_socp`m'" "th_comp`m'" "th_othp`m'"    { 
/*capture:*/ svmat _`coef' 		// store average markdown for this bootstrap iteration in matrix
} 
}
 forvalues m = 1/2 {
foreach coef in    "th_union`m'"    { 
/*capture:*/ svmat _`coef' 		// store average markdown for this bootstrap iteration in matrix
} 
 
}
 
keep in 1/`B'
list _bl* _bm* _bk*

* Reshape and summarize kappa by year. Extract SD per year (or 5%, 95% CIs per year)? Plot this.
   
save ./data/temp/bootstrapsap, replace
use ./data/temp/bootstrapsap, clear

list m_md_nu5, d

sum m_kappa*, d
sum _kappa*, d

sum _rho1*
 

sum _bl_bb2*

  sum _bl_bb11, d

  sum _bm_nu5, d
  sum _bm_bb11, d

  
 

** Standard errors and 5%-95% C.I.s

foreach var of varlist *_md*  *_mu*  _th*   _b*  _scale*  *_kappa* *_rho* {
egen se`var' = sd(`var')
egen p5`var' = pctile(`var'),p(5)
egen p95`var' = pctile(`var'),p(95)
}
  

 
  
 
** Production function by time block

 
 
foreach mod in "bl" "bm" "bk" { 
rename (se_`mod'_blockbis1 se_`mod'_blockbis2  )(b1 b2  )
estpost su b1 b2   
est store `mod'_timev1_se   
drop b1 b2  
}
 
foreach mod in "bl" "bt" "blt" { 
rename (se_`mod'_bb5 )(b)
estpost su b
est store `mod'_timev3_se   
drop b
}
     

* IV comparisons 


foreach var in "l" "k" "m" {
forvalues n = 1/3 { 
gen b_iv`n' = se_b`var'_iv`n' 
}
estpost su b_iv1 b_iv2 b_iv3 
drop b_iv* 
est store b`var'_iv_se 
}

 * returns to scale 
 
forvalues n = 0(5)15 { 
 foreach mod in "bl" "bm" "bk" { 
rename (se_`mod'_nu`n'   )(`mod'_nu  )
}
 foreach mod in "md" "mu" { 
rename (sem_`mod'_nu`n'   )(`mod'_nu  )
}
estpost su bl_nu bm_nu bk_nu md_nu mu_nu
est store pf_nu`n'_se 
drop bl_nu bm_nu bk_nu md_nu mu_nu 
}
	 
 
*  translog and interaction term model
 
foreach mod in "bl" "bm" "bk" "blm" "bmk" "blk" "bl2" "bm2" "bk2"  { 
rename (se_`mod'_bb3   )(`mod'_bb3  )
}
estpost su bl_bb3  bm_bb3 bk_bb3 blm_bb3 blk_bb3 bmk_bb3 bl2_bb3 bm2_bb3 bk2_bb3
est store pf3_bb_se
 
* CS-based coefficient SEs
 
foreach var in "bl" "bm" "bk" "mu" { 
rename (se_`var'_cs_col se_`var'_cs_noc) (col noc) 
estpost su col noc  
est store `var'_cs_se  
drop col noc
}

* Markup response to cartel

forvalues m = 1/2 {
rename (se_mucar`m') (mucar) 
estpost su mucar
est store mucar`m'_se  
drop mucar
}

* Rev pf

forvalues m = 6/11 {
rename (se_bl_bb`m'1 se_bm_bb`m'1 se_bk_bb`m'1 se_rho_bb`m'1) (bl bm bk rho) 
estpost su bl bm bk rho
est store pf_bb`m'_se  
drop bl bm bk rho
}


forvalues m = 1/1 {
rename (se_bl_bb`m'1 se_bm_bb`m'1 se_bk_bb`m'1 se_rho_bb`m'1) (bl bm bk rho) 
estpost su bl bm bk rho
est store pf_bb`m'_se  
drop bl bm bk rho
}


forvalues m = 12/12 {
rename ( se_bl_bb`m' se_bm_bb`m' se_bk_bb`m' se_rho1_bb`m'  se_rho2_bb`m') (bl bm bk rho1 rho2) 
estpost su bl bm bk rho1 rho2
est store pf_bb`m'_se  
drop bl bm bk rho1 rho2
}
 
 
forvalues m = 1/12 {
 rename (se_mu_bb`m'1 se_md_bb`m'1  ) (mu md) 
 estpost su mu md 
 est store mu_bb`m'_se  
 drop mu md 
}

rename se_bnu bnu 
estpost su bnu
est store bnu_se 
drop bnu 

  
** Store results to table

* interaction effects
   
local m = 2
foreach mod in "ols" "bb" {
rename (se_bl_`mod'`m' se_bk_`mod'`m'  se_bm_`mod'`m'   se_blk_`mod'`m' 	 se_bmk_`mod'`m'   se_thl_`mod'`m' se_thm_`mod'`m'	se_scale_`mod'`m'   ) (bl bk bm  blk bmk thl thm scale  )
estpost su bl bm bk   blk bmk
est store pf_`mod'`m'_se
estpost su thl thm scale  
est store th_`mod'`m'_se
drop bl bm bk blk  thm bmk thl scale
}
  
  
